!dk mai
!
      program celest3d
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
!
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv
      integer , dimension(4) :: lioinput, lioprint
      integer :: lioplot, lbox, n, nsm, ll, nn, iter
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5
      real(double), dimension(8) :: wght
      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
         ucm, vcm, wcm, cm
      real(double), dimension(3,20) :: wsin, wcos
      real(double), dimension(3,12,20) :: tsi, rot
      real(double) :: sok, sko, cok, cko, celplt, wk, absbc, rvvol, where, &
         transpos, trl, tmin, told, xx, t2old, error, relax, tiny, dumdt, phibc&
         , dumscal, ws, errormaxwell, t3
      logical :: expnd, nomore, firstime,divergence_cleaning
      real(double) :: dnorm,ez_avg !,eenbefore,eenafter
!-----------------------------------------------
!
!
!------------------------------------------------------------------
!               main program
!------------------------------------------------------------------
!
!     make calls to timing routines
!
!
      divergence_cleaning=.true.
!
      open(unit=5, file='input', position='asis')
!ll   1. open i/o devices:
!        ----------------
!      call link ("unit22 = (input   ,open),
!     .            unit9  = (out1 ,create,text),
!     .            unit6  = (out2 ,create,text),
!     .            unit99 = (film    ,create,text),
!     .            unit15 = (pldata,create,len=10b),
!     .            unit16 = (box     ,create,len=10b),
!     .            //")
!
      lioinput(1) = 22
      lioinput(2) = 0
      lioinput(3) = 0
      lioinput(4) = 0
      lioprint(1) = 9
      lioprint(2) = 6
      lioprint(3) = 99
      lioprint(4) = 0
      lioplot = 15
      lbox = 16
!
!
!l    krd    ...device number for read input cards
!l    kpt    ...device number for print
!l    kpr    ...device number for print (switched)
!l    kpl    ...device number for storing data on disk
!l    kfm    ...device number for film/microfiche
!l    ktpin  ...device number for binary tape read
!l    ktpout ...device number for binary tape write
      krd = 5
      kpt = 6
      kpl = 15
      kfm = 59
      ktpin = 7
      ktpout = 7
!
!ll   2. initialize run:
!        --------------
!
!     initialize control counters
!
      sok = 0
      sko = 0
      cok = 0
      cko = 0
!

 	  pixxsp2 = 0.0
          pixysp2 = 0.0
          pixzsp2 = 0.0
          piyysp2 = 0.0
          piyzsp2 = 0.0
          pizzsp2 = 0.0
 	  pixxsp1 = 0.0
          pixysp1 = 0.0
          pixzsp1 = 0.0
          piyysp1 = 0.0
          piyzsp1 = 0.0
          pizzsp1 = 0.0
!
          densp1=0d0
		  uxsp1=0d0
		  uysp1=0d0
		  uzsp1=0d0
!
          densp2=0d0
		  uxsp2=0d0
		  uysp2=0d0
		  uzsp2=0d0
!
          exavp = 0.0
          eyavp = 0.0
          ezavp = 0.0
!
          avg=0d0
!
          potavp = 0.0
!
!     open graphics library
!
      call gplot ('u', celplt, 6)
      call libdisp
      write (0, *) ' celest3d: starting initiv'

      call initiv

      if (restart) then
         call restartin
         call plotinit (ibp2, jbp2, kbp2, nrg)
         ncyc1 = 1
         ncyc2 = ncyc_restart
         write (*, *) 'ncyc1,2', ncyc1, ncyc2
         call date_and_time(ddate,hour)
         write(*,*)'date (ccyymmdd) = ',ddate
         write(*,*)'time (hhmmss.sss) = ',hour
         go to 999
      endif
      call plotinit (ibp2, jbp2, kbp2, nrg)
!
!     *******************************************************
!
!     calculate the gradient of the magnetic field strength only if
!     the drift approximation is used
!
!     **********************************************************
      if (anydrift) then
!
         do k = 1, kbp2
            do j = 1, jbp2
               i = 1
               if (ibp2 > 0) then
!
                  absb(1+(j-1)*jwid+(k-1)*kwid:(ibp2-1)*iwid+1+(j-1)*jwid+(k-1)&
                     *kwid:iwid) = 0.0
!
                  i = ibp2 + 1
                  ijk = (ibp2 - 1)*iwid + 1 + (j - 1)*jwid + (k - 1)*kwid
               endif
            end do
         end do
!
         write (0, *) 'celest3d: starting absolute'
         call absolute (ncells, ijkcell, bxn, byn, bzn, absb)
!c
         wk = 0.
         absbc = 0.
         call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkr, wkr, wkf, wke, &
            absbc, absb, periodicx)
!
!
         call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
            c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, absb, gradxb, gradyb, gradzb)
!
         call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradxb, gradyb, gradzb)
!
!
         do n = 1, nvtx
!
            ijk = ijkvtx(n)
!
            rvvol = 1./vvol(ijk)
!
            gradxb(ijk) = gradxb(ijk)*rvvol
            gradyb(ijk) = gradyb(ijk)*rvvol
            gradzb(ijk) = gradzb(ijk)*rvvol
!
         end do
!
      endif
!
!
!     ************************************************
!
!      initialize the particles
!
!     ************************************************
!
      write (0, *) ' celest3d: starting parset'
      call current
      call parset
      if (nrg_obj > 0) call parset_obj
!
      where = 0.
!      call counter(ncells,ijkcell,iphead,link,number,npart,where)
!
!     **************************************************
!
!     interpolate particle quantities to the grid
!
!     **************************************************
!
      write (0, *) ' celest3d: starting parcelc'
      call parcelc
!
!
      transpos = 1.
!
      call parcelv (ncells, ijkcell, nvtx, ijkvtx, nsp, iwid, jwid, kwid, ibar&
         , jbar, kbar, dt, ijkctmp, anydrift, eandm, itdim, iphead, iphd2, link&
         , transpos, qom, ico, qpar, pxi, peta, pzta, up, vp, wp, pmu, bxv, byv&
         , bzv, wate, uptilde, vptilde, wptilde, mask, pixx, pixy, pixz, piyy, &
         piyz, pizz, qdnv, number, jx0, jy0, jz0, jmag, jxtilde, jytilde, &
         jztilde, jxs, jys, jzs,  pixysp2, piyzsp2, pixzsp2, pixxsp2, piyysp2, pizzsp2, &
         pixysp1, piyzsp1, pixzsp1, pixxsp1, piyysp1, pizzsp1, uxsp1,uysp1,uzsp1,densp1,uxsp2,uysp2,uzsp2,densp2, 0d0)
!
!
!     ******************************************************
!
!     interpolate object quantities to the grid
!
!     *****************************************************
!
      if (nrg_obj > 0) then
!
         write (*, *) 'celest3d: starting parcelv_obj'
!
         call parcelv_obj (ncells, ijkcell, nvtx, ijkvtx, nrg_obj, iwid, jwid, &
            kwid, ijkctmp, itdim, iphead_obj, iphd2, link_obj, ico_obj, &
            massp_obj, pxi_obj, peta_obj, pzta_obj, chip_obj, wate, chi, mv)
!
         write (*, *) 'celest3d: starting parcelc_obj'
!
         call parcelc_obj
      endif
!
!
      write (*, *) 'celest3d: starting vinit'
!
      call vinit
!
      write (*, *) 'celest3d: starting budget'
!
      call budget
!
!     ******************************************************************
!
!     generate an adaptive grid
!
!     *****************************************************************
!
      if (adaptg) then
!
         call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
            c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, chic, gradxchi, gradychi, gradzchi)
!
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, gradxchi, gradychi, gradzchi, periodicx)
!
         call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradxchi, gradychi, &
            gradzchi)
!
!
!     set gradient at k=kbp2 equal to gradient at k=kbp1
!
         do j = 1, jbp2
            do i = 1, ibp2
!
               ijk = (kbp2 - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1
!
               gradxchi(ijk-kwid) = gradxchi(ijk)
               gradychi(ijk-kwid) = gradychi(ijk)
               gradzchi(ijk-kwid) = gradzchi(ijk)
!
            end do
         end do
!
         do n = 1, nvtx
!
            ijk = ijkvtx(n)
!
            chic(ijk) = sqrt(gradxchi(ijk)**2+gradychi(ijk)**2+gradzchi(ijk)**2&
               )/vvol(ijk) + 1.
!
         end do
!
         write (*, *) 'celest3d: calling smooth'
!
         nsm = 10
!
         call smooth (2, ibp2, 2, jbp2, 2, kbp2, iwid, jwid, kwid, nsm, chic, &
            gradxchi, periodicx)
!
         write (*, *) 'celest3d: starting meshgen'
!
         call meshgen (ibp2, jbp2, kbp2, x, y, z, wgrid, chic, vol, dt, taug, &
            numgrid, wate(1,7), wate(1,8), wate(1,9), wate(1,10), wate(1,11), &
            wate(1,12), wate(1,13), wate(1,14), wate(1,15), lama, lamb, lamc, &
            wate(1,1), wate(1,2), wate(1,3), delta, psi, tsix, tsiy, tsiz, etax&
            , etay, etaz, nux, nuy, nuz, wate(1,4), wate(1,5), wate(1,6), iwid&
            , jwid, kwid, ijktmp2, istep, iota, periodic, toroid, strait, rmaj&
            , rwall, btorus, q0, cdlt, sdlt, dz)
!
         write (*, *) 'celest3d: starting celindex'
!
         call celindex (2, ibp1, 2, jbp1, 2, kbp1, iwid, jwid, kwid, ijkcell, &
            ijkctmp, ncells, ijkvtx, nvtxkm, nvtx)
!
         write (*, *) 'celest3d: calling geom'
!
         call geom (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibar, jbar&
            , kbar, nvtxkm, cartesian, cdlt, sdlt, dz, x, y, z, c1x, c2x, c3x, &
            c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
            c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, vvol_half, &
            ijktmp2, a11, a12, a13)
!
!
         write (*, *) 'celest3d: starting metric'
!
         call metric (ncells, ijkcell, iwid, jwid, kwid, x, y, z, tsix, tsiy, &
            tsiz, etax, etay, etaz, nux, nuy, nuz, sie)
!
!     relocate particles on the new mesh
!
         write (*, *) 'celest3d: starting parlocat_all'
!
         call parlocat_all
         write (*, *) 'celest3d: starting parlocat_all_obj'
!
!
         call parlocat_all_obj
!
         write (*, *) 'celest3d: starting parcelv_obj'
!
         call parcelv_obj (ncells, ijkcell, nvtx, ijkvtx, nrg_obj, iwid, jwid, &
            kwid, ijkctmp, itdim, iphead_obj, iphd2, link_obj, ico_obj, &
            massp_obj, pxi_obj, peta_obj, pzta_obj, chip_obj, wate, chi, mv)
!c
         call parcelc
!c
!     ADDED BY G. LAPENTA
         call parcelv (ncells, ijkcell, nvtx, ijkvtx, nsp, iwid, jwid, kwid, &
            ibar, jbar, kbar, dt, ijkctmp, anydrift, eandm, itdim, iphead, &
            iphd2, link, transpos, qom, ico, qpar, pxi, peta, pzta, up, vp, wp&
            , pmu, bxv, byv, bzv, wate, uptilde, vptilde, wptilde, mask, pixx, &
            pixy, pixz, piyy, piyz, pizz, qdnv, number, jx0, jy0, jz0, jmag, &
            jxtilde, jytilde, jztilde, jxs, jys, jzs,  &
            pixysp2, piyzsp2, pixzsp2, pixxsp2, piyysp2, pizzsp2, &
            pixysp1, piyzsp1, pixzsp1, pixxsp1, piyysp1, pizzsp1, uxsp1,uysp1,uzsp1,densp1,uxsp2,uysp2,uzsp2,densp2,0d0)
!
         call vinit
         call output
         call budget
!
      endif
!
!l    (write data file on disk)
!      call pltset (0)
!
!      call second (stim)
      write(*,*)'second-tlm',tlm
      call getjtl (itlm)
           write(*,*)'getjtl',itlm
      itlwd = 0
      t1 = stim
      t2 = t1
      trl = t1
      tmin = (0.02 + 0.005)*fibar*fjbar*fkbar
      tlm = itlm - tmin
      if (tlm <= 0) then
!        (time limit too small for problem)
         write (kpt, 9020) tmin
         write (kfm, 9020) tmin
         call exit (1)
!
!
!ll   3. time-step: integrate equations in time:
!        --------------------------------------
      endif
!l       (compute normal to outer surface)
!         call normal
!
  999 continue
      if (ncyc2 > 0) then
!
!
!
!l    >>>>>>>>>>>>>>> ncyc-loop <<<<<<<<<<<<<<<
!
         do ncyc = ncyc1, ncyc2
         !
            nh = mod(ncyc,nhst)
            nh = max(1,min(nhst,nh))
            t = t + dt
            write (*, *) 'TIME', ncyc, t, ncyc2
            told = t2
            !call second (t2)
            xx = (t2 - told)*rbjbkb
            if (ncyc == ncyc1) then
               costcc(nh) = 0.
            else
               costcc(nh) = (t2 - t2old)/(ibar*jbar*kbar)
            endif
            t2old = t2
!
!     ****************************************************
!
!     prepare for the next computation cycle
!
!     ****************************************************
!

            if (.not.testpar .and. (ncyc>1 .or. restart)) then
!
               call vinit
!
            endif

!
            call budget
!
            if (lpr /= 0) then
!
               if (ncyc - 1==0 .or. mod(ncyc,50)-1==0) pflag = 1.
               if (pflag > 0.) write (kpt, 9000)
               write (kpt, 9010) ncyc, t, tmass, tmomx, tmomy, tmomz, racc, &
                  zacc, tie, tke, tbe, tle, numitv, numit, dvbmax, xx
               pflag = -1.
!
!
!      if(mod(ncyc,nfreq).eq.0) call output
               call output
!
            endif
!     <<<<<<<<<<<<<<<<<<<<<< Solve Maxwell's Equations >>>>>>>>>>>>>>>>>>>>>>>>>
!
!
!     zero working arrays
!

            wate(:itdim,:5) = 0.0
            phipot(:itdim) = 0.0
!
            if (fields) then
               error = 1.E-3
!               error=1.e-2
               relax = 1.
               tiny = 0.E-2

!

	if(divergence_cleaning) then

!     calculate the source for Poisson's equation
!     the call computes a correction to E0
!
               call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x&
                  , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y&
                  , c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, ex0, ey0, ez0&
                  , divpix)
!
	dnorm=0.
               do n = 1, ncells
                  ijk = ijkcell(n)
                  wate(ijk,10) = fourpi*qdnc0(ijk) - divpix(ijk)
		  dnorm=dnorm+(fourpi*qdnc0(ijk) - divpix(ijk))**2
               end do
	       write(*,*)'prima=',dnorm
		   !eenbefore=sum(ex0(ijkvtx(1:nvtx))**2)+sum(ey0(ijkvtx(1:nvtx))**2)+sum(ez0(ijkvtx(1:nvtx))**2)
!
!
!     with these parameters, poisson solves poisson's equation
!
               dumdt = 0.0
               dumscal = 1.
               wk = 0.
               phibc = 0.
!
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,1), wate(1,2)&
                  , wate(1,7), wate(1,8), itmax, error, wate(1,10), dumdt, &
                  dumscal, wk, phibc, ex, ey, ez, wate(1,4), wate(1,5), wate(1,&
                  6), phipot, wate(1,9), wkl, wkr, wkf, wke, periodicx)
!
!
               call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, &
                  c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
                  c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, phipot, gradxb, &
                  gradyb, gradzb)
!
!
               call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradxb, gradyb, &
                  gradzb) !....nothing is done with this call
!
!dir$ ivdep
               !ez_avg=0d0
               do n = 1, nvtx
                  ijk = ijkvtx(n)
                  rvvol = 1./vvol(ijk)
                  ex0(ijk) = ex0(ijk) - gradxb(ijk)*rvvol
                  ey0(ijk) = ey0(ijk) - gradyb(ijk)*rvvol
                  ez0(ijk) = ez0(ijk) - gradzb(ijk)*rvvol
                  !ez_avg=ez_avg+ez0(ijk)
               end do
               !ez0=ez0-ez_avg/dble(nvtx)

         call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x&
                  , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y&
                  , c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, ex0, ey0, ez0&
                  , divpix)
!
			   dnorm=0.
               do n = 1, ncells
                  ijk = ijkcell(n)
                  wate(ijk,10) = fourpi*qdnc0(ijk) - divpix(ijk)
				  dnorm=dnorm+(fourpi*qdnc0(ijk) - divpix(ijk))**2
               end do
	       write(*,*)'dopo=',dnorm
		   !eenafter=sum(ex0(ijkvtx(1:nvtx))**2)+sum(ey0(ijkvtx(1:nvtx))**2)+sum(ez0(ijkvtx(1:nvtx))**2)
		   !write(*,*)'Jerrys puzzle ',eenbefore,eenafter
		   !write(101,*)'before=',eenbefore,'    after=',eenafter

	end if
!
!	Applying BC to the new E0
!
               call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, &
                  sdlt, cdlt, ex0, ey0, ez0, periodicx)
!
!	Smmothing the new E0
!
               phibc = 0.
               call vtoctov (nvtx, ijkvtx, iwid, jwid, kwid, ibp2, jbp2, kbp2, &
                  nsmooth, sdlt, cdlt, wkl, wkr, wkf, wke, phibc, ex0, ey0, ez0, &
                  a11, a12, a13, periodicx)

!
!
!    Electrostatic  potential
!
!
!
               do n = 1, ncells
                  ijk = ijkcell(n)
                  wate(ijk,10) = fourpi*qdnc0(ijk)
               end do
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,1), wate(1,2)&
                  , wate(1,7), wate(1,8), itmax, error, wate(1,10), dumdt, &
                  dumscal, wk, phibc, ex, ey, ez, wate(1,4), wate(1,5), wate(1,&
                  6), phipot, wate(1,9), wkl, wkr, wkf, wke, periodicx)
!
               potavp(:ibp2*jbp2*kbp2) = potavp(:ibp2*jbp2*kbp2)*.99d0&
                    + phipot(:ibp2*jbp2*kbp2)*0.01d0
!
!
               if (eandm) then
!
!     calculate the source for Faraday's law
!     store in jxtildet,jytildet,jztildet
!
!
                  ws = fourpi/clite
                  dtc = clite*cntr*dt
!
!
!
!
                     call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x&
                        , c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y&
                        , c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, &
                        qdnc, jxtildet, jytildet, jztildet)
!
!
!
!
                  call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, &
                     ijkvtx, vol, wate(1,28))
!
!     wk=1.
!     phibc=0.
!      call vtoctov(nvtx,ijkvtx,iwid,jwid,kwid,ibp2,jbp2,kbp2,
!     &     nsmooth,sdlt,cdlt,wkl,wkr,wkf,wke,phibc,
!     &jxtilde,jytilde,jztilde,a11,a12,a13,periodicx)

!     &jx0,jy0,jz0,a11,a12,a13)
!

                     do n = 1, nvtx
                        ijk = ijkvtx(n)
!
                        rvvol = 1./vvol(ijk)
                        ws = fourpi/clite*vvol(ijk)
!
                        jxtildet(ijk) = (ex0(ijk)-fourpi*jxtildet(ijk)*rvvol*&
                           dtc*dtc+dtc*(jx0(ijk)-ws*jxtilde(ijk))*rvvol)*wate(&
                           ijk,28)
!
                        jytildet(ijk) = (ey0(ijk)-fourpi*jytildet(ijk)*rvvol*&
                           dtc*dtc+dtc*(jy0(ijk)-ws*jytilde(ijk))*rvvol)*wate(&
                           ijk,28)
!
                        jztildet(ijk) = (ez0(ijk)-fourpi*jztildet(ijk)*rvvol*&
                           dtc*dtc+dtc*(jz0(ijk)-ws*jztilde(ijk))*rvvol)*wate(&
                           ijk,28)
!
                     end do
!
!write(*,*)'ex0',sum(ex0**2),sum(ey0**2),sum(ez0**2)
!write(*,*)'ex0',sum(jxtilde**2),sum(jytilde**2),sum(jztilde**2)
!jxtildet=jx0
!jytildet=jy0
!jztildet=jz0
!
                     ex(:ibp2*jbp2*kbp2) = ex0(:ibp2*jbp2*kbp2)
                     ey(:ibp2*jbp2*kbp2) = ey0(:ibp2*jbp2*kbp2)
                     ez(:ibp2*jbp2*kbp2) = ez0(:ibp2*jbp2*kbp2)
!
                     errormaxwell = error
!
                     call maxwell (nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, &
                        precon, tiny, ncells, ijkcell, itdim, nrg, ibp1, jbp1, &
                        kbp1, wate(1,28), wate(1,29), cdlt, sdlt, strait, dx, dy, dz, &
                        clite, cntr, dt, t, fourpi, itmax, errormaxwell, &
                        periodicx, bx_ini, c1x, &
                        c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, &
                        c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
                        c8z, vvol, vol, x,y,z,qom, qdnv, qdnc, curlx, curly, curlz, wate(1,&
                        1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                        6), wate(1,7), wate(1,8), wate(1,9), wate(1,10), wate(1&
                        ,11), wate(1,12), wate(1,13), wate(1,14), wate(1,15), &
                        wate(1,16), wate(1,17), wate(1,18), wate(1,19), wate(1,&
                        20), wate(1,21), a11, a12, a13, a21, a22, a23, a31, a32&
                        , a33, iphd2, wate(1,22), wate(1,23), wate(1,24), wate(&
                        1,25), wate(1,26), wate(1,27), jxtilde, jytilde, &
                        jztilde, jx0, jy0, jz0, jxtildet, jytildet, jztildet, &
                        pixx, pixy, pixz, piyy, piyz, pizz, divpix, divpiy, &
                        divpiz, colx, coly, colz, udrft, vdrft, wdrft, tsix, &
                        tsiy, tsiz, bxv, byv, bzv, ex0, ey0, ez0, ex, ey, ez, &
                        iter,ex_ext,ey_ext,ez_ext,shock,susxx,susxy&
                        ,susxz,suszx,suszy,suszz)
!
                  citer_maxwell(nh) = iter

         call curlv (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, bxn, byn, bzn, jx0, jy0, jz0)
!
            endif
!
            phibc = 0.
            wk = 1.
!
!
               call vtoctov (nvtx, ijkvtx, iwid, jwid, kwid, ibp2, jbp2, kbp2, &
                  0, sdlt, cdlt, wkl, wkr, wkf, wke, phibc, ex, ey, ez, a11, &
                  a12, a13, periodicx)
!
!
!
!     **************************************************
!
!	add an external electric field or other force
!	mapped into an effective electric field
!
!     **************************************************
!
!	call external_field

	endif
!
!
!
!     **************************************************
!
!     solve the particle equations of motion
!
!     **************************************************
!
!
            if (testpar) then
!
               ncellsp = 1
               i = int(pxi(nphist+1))
               j = int(peta(nphist+1))
               k = int(pzta(nphist+1))
!
               ijkctmp(1) = (k - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1
!
            else
!
               ncellsp = ncells
!
            endif
!
            if (anydrift) then
!
               where = 1.
!      call counter(ncells,ijkcell,iphead,link,number,npart,where)
               call parmovgc
               where = 1.1
!      call counter(ncells,ijkcell,iphead,link,number,npart,where)
            endif
!
            if (explicit) then
!
!      call parmove
               if(relativity) then
               call parmovi_rel
               else
               call parmovi
               end if
!
            else
!

!
               if(relativity) then
               call parmovi_rel
               else
               call parmovi
               end if

!
!
!     particle control
!
     call control(sok,sko,cok,cko)
     write(*,*)'splitting: fail=',sko,'    succ=',sok
     write(*,*)'coalescen: fail=',cko,'    succ=',cok


            endif
!
!
            call ParInject
!
            if (testpar) call celindex (2, ibp1, 2, jbp1, 2, kbp1, iwid, jwid, &
               kwid, ijkcell, ijkctmp, ncells, ijkvtx, nvtxkm, nvtx)
!
!     **************************************************
!
!     interpolate particle quantities to the grid
!
!     **************************************************
!
            if (.not.testpar) then
               where = 5.
!      call counter(ncells,ijkcell,iphead,link,number,npart,where)
               call parcelc
!
!
               transpos = 1.
!
               call parcelv (ncells, ijkcell, nvtx, ijkvtx, nsp, iwid, jwid, &
                  kwid, ibar, jbar, kbar, dt, ijkctmp, anydrift, eandm, itdim, &
                  iphead, iphd2, link, transpos, qom, ico, qpar, pxi, peta, &
                  pzta, up, vp, wp, pmu, bxv, byv, bzv, wate, uptilde, vptilde&
                  , wptilde, mask, pixx, pixy, pixz, piyy, piyz, pizz, qdnv, &
                  number, jx0, jy0, jz0, jmag, jxtilde, jytilde, jztilde, jxs, &
                  jys, jzs,  pixysp2, piyzsp2, pixzsp2, pixxsp2, piyysp2, pizzsp2, &
                  pixysp1, piyzsp1, pixzsp1, pixxsp1, piyysp1, pizzsp1, &
				  uxsp1,uysp1,uzsp1,densp1,uxsp2,uysp2,uzsp2,densp2,avg)
!
!
               if (nrg_obj > 0) then
!
                  call parcelv_obj (ncells, ijkcell, nvtx, ijkvtx, nrg_obj, &
                     iwid, jwid, kwid, ijkctmp, itdim, iphead_obj, iphd2, &
                     link_obj, ico_obj, massp_obj, pxi_obj, peta_obj, pzta_obj&
                     , chip_obj, wate, chi, mv)
!
               endif
!
!
!
            endif
!
!

            !call second (t3)
            t3=0d0
!            This line below used to stop the simulation when a certain total simulation
!            time was reached (actual cpu wallclock)
!            if (t>=twfin .or. t3-stim>=tlm) go to 1010

!
            if (t2 - t1 < dcttd) then
!
               if (t + 1.E-10 < pttd) cycle
!
               pttd = pttd + dpttd
!
            endif
!l       (write all necessary data on tape (dump))
            t1 = t2
!
         end do
!
      endif
!
!l    >>>>>>>>>>>>>>> end of ncyc-loop <<<<<<<<<<<<<<<
!
      ncyc = ncyc2
!
!ll   4. terminate run:
!        -------------
      call output
!
 1010 continue
  212 format(i5,9(1p,e11.3))
!
!
      call finito
!
!
      call plotclose
!
      call restartout
!
      write (*, *) 'Im stopping from celest3d after ncyc=', ncyc
      stop
!
 9000 format('  cyc',8x,'t',6x,'mass',4x,'x-mom',3x,'y-mom',3x,'z-mom',6x,&
         'racc',6x,'zacc',5x,'dei',5x,'dek',5x,'deb',5x,'etotl',' itv',' it',2x&
         ,'dvbmax',4x,2x,'grinds')
 9010 format(1x,i4,1p,e9.2,1e11.3,1x,3e8.1,2e10.3,3e8.1,e10.3,i4,i3,e8.1,4x,e&
         8.1)
 9020 format('>>> error: time limit must exceed ',1p,e9.2,'seconds <<<')
      stop
!
      end program celest3d

